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We present molecular dynamics simulations of a binary Lennard- Jones mixture at temperatures 
below the kinetic glass transition. The "mobility" of a particle is characterized by the amplitude of its 
fluctuation around its average position. The 5% particles with the largest /smallest mean amplitude 
are thus defined as the relatively most mobile/immobile particles. We investigate for these 5% particles 
their spatial distribution and find them to be distributed very heterogeneously in that mobile as well as 
immobile particles form clusters. The reason for this dynamic heterogeneity is traced back to the fact 
that mobile/immobile particles are surrounded by fewer/more neighbors which form an effectively 
wider/narrower cage. The dependence of our results on the length of the simulation run indicates 
that individual particles have a characteristic mobility time scale, which can be approximated via the 
non- Gaussian parameter. 
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I. INTRODUCTION 

Although glasses have already been studied for a long 
time, their complete understanding is still an open prob- 
lem due to the complex behavior of their static and dy- 
namic quantities jy. We focus here on their dynamics 
which has been found to relax non-exponentially in the 
supercooled liquid and shows strong history dependence 
below the glass transition. 

The question arises whether this behavior is due to spa- 
tially homogeneous non-exponential dynamics or spatially 
heterogeneous dynamics (for review articles see 
Since a glass is an amorphous solid and therefore not all 
atoms arc structurally equivalent, one should expect that 
also their dynamics differs, i.e. that the system has a het- 
erogeneous dynamics. Recently the answer to the ques- 
tion of dynamic heterogeneity has been addressed both 
by means of experiments ]5|-|7j and computer simulations 
of two-dimensional p|-[T()[ and three-dimensional systems 
pl|-[l^]. Here we study a binary Lennard- Jones mixture 
in three dimensions which has been investigated before 
extensively [ [TBI and which shows clear dynamic hetero- 
geneity above the glass transition. Similar dynamics has 
been found experimentally with confocal microscopy of a 
colloidal suspension in the supercooled fluid and in the 
glass [^,0. Whereas most simulations have been done at 
relatively high temperatures above the calorimetric glass 
transition, and the experiments of atomic systems were 
performed near the glass transition, we simulate, in this 
paper, below the glass transition (which has so far only 
been done experimentally by Weeks et al. and in simu- 
lations by Oligschleger et al. |l6|]). We find here dynamic 
heterogeneity via simulations of the same binary Lennard- 
Jones system as but below the glass transition. In con- 
trast to the previous simulations we have the picture of a 
solid in mind, instead of coming from the liquid. We use 



the "localization length" of the work of reference |17| to 
define the mobility of a particle as the mean fluctuation 
around its average position. To address the question of 
what allows or inhibits a particle to be mobile, we study 
the surrounding of these particles. Using different lengths 
of simulation runs we also learn about the time scale over 
which mobile and immobile particles sustain their charac- 
ter. 

We review in section |0] the model used and give de- 
tails of the simulation. In section [II we present the mean 



square displacement and the mean fluctuations of a parti- 
cle around its average position and define what we mean 
by mobile and immobile particles. We then study their 
spatial distribution (sec. IV), surrounding (se c. |v| ) and 
time scale (sec. VI ) and conclude with section |VII| . 



II. SIMULATION DETAILS 

We study a binary Lennard- Jones (LJ) mixture of 800 
A and 200 B particles. Both A and B particles have the 
same mass. The interaction between two particles a and 
p (a, /3 g {A,B}) is 



(1) 



where eaa = L0, cab = L5, cbb = 0.5, ctaa = 1-0, (t A b = 
0.8 and ctbb = 0.88. We truncate and shift the potential 
at r = 2.5a a f3. From previous investigations JLsfl it is 
known that this system is not prone to crystallization and 
demixing. In the following we will use reduced units where 
the unit of length is <taa, the unit of energy is 6aa and 
the unit of time is yj ma\ A / (48eaa) ■ 

We carry out molecular dynamics (MD) simulations us- 
ing the velocity Verlet algorithm with a time step of 0.02. 
The volume is kept constant at V = 9.4 3 = 831 and we 



1 



use periodic boundary conditions. We are interested in 
the dynamics of the system below the glass transition. 
Since recent simulations Jl3] showed that for present day 
computer simulations the system falls out of equilibrium 
around T w 0.44, we run (NVE)-simulations at tempera- 
tures T = 0.15/0.2/0.25/0.3/0.35/0.38/0.4/0.41/0.42 and 
0.43. To do so we start with a well equilibrated config- 
uration at T = 0.466. After an instantaneous quench to 
T = 0.15 we first run a (NVT)-simulation Q for 10 5 MD 
steps to let the system anneal, and then run the produc- 
tion run with a (NVE) simulation for also 10 5 MD steps. 
We then increase the temperature to T = 0.2 and then 
again run a (NVT)-simulation followed by a production 
run each with 10 5 MD steps; then increase the tempera- 
ture to T — 0.25 and so forth. In this paper we refer to 
the so obtained production runs as "short runs," of which 
some preliminary results have been published elsewhere 
. We present here mainly results of the so-called "long 
runs" for which we use the configurations at the end of the 
equilibration period of the short runs but the production 
runs are for 5 • 10 6 MD steps. To improve the statistics we 
run 10 independent configurations for both long and short 
runs and for each temperature. 

As it has been demonstrated in earlier work |pl| , the 
structural properties of glasses studied by molecular dy- 
namics simulation do depend on the preparation history 
quite distinctly. Since we study the system out of equi- 
librium and at finite temperature, the resulting config- 
urations show some "aging phenomena" during the time 
intervals used for the production of the present data. How- 
ever, for those temperatures where even the mean square 
displacements of the 5% fastest particles are still small 
in comparison to typical interparticle distances over the 
whole time span of the averaging, this change of the glass 
structure due to aging should have a relatively small ef- 
fect on our data. As will be seen below, this is the case for 
T < 0.35, while for somewhat higher temperatures aging 
effects can be expected to become important. The reason 
for this is that r, the typical relaxation time of the system 
in equilibrium, is at T — 0.446 around 800, 000 time units 
(= 4- 10 7 MD steps) |§], thus only a factor 10 longer than 
the long runs in the present work. Since in equilibrium the 
a— relaxation time corresponds to the typical time scale on 
which a substantial fraction of the particles have moved 
(40-70%), it can be expected that quite a few particles will 
show relaxation even on time scales significantly shorter 
than r. Similar relaxation processes are also expected in 
the out of equilibrium situation, i.e. in the glass and thus 
we do indeed expect aging effects at temperatures slightly 
below the (kinetic) glass transition. We shall comment on 
this problem that there is some aging of the glass structure 
occurring when appropriate. 



III. MOBILE AND IMMOBILE PARTICLES 

Similar to previous work on dynamic heterogeneities 
[^ ]l4| , |ll| , we study the dynamics of the system by ob- 
serving the fastest (mobile) and the slowest (immobile) 
particles. Since the focus of this work is, however, on the 
dynamics of the glass below the glass transition, our def- 
inition of mobile and immobile particles is different. We 
have in mind the picture of a harmonic solid for which the 
vibrational amplitude carries essential information about 
the local dynamics. We therefore characterize the mobility 
of each particle i by 



d = 



— ,2 



(2) 



where the bar denotes an average over a certain time inter- 
val. We call the 5% A particles and separately the 5% B 
particles with the largest /smallest d\ the mobile/immobile 
particles. With "mobile" we intend to indicate that these 
particles are relative to all other particles more mobile, 
they are however in the results, presented here, in most 
cases still bound to their site. All results presented be- 
low are qualitatively the same if the 5% are replaced by 
10% particles and are therefore independent of the spe- 
cific percentage used. The results depend however on the 
length of the time average, as will be discussed in detail in 
section V] . We use in this paper the non-Gaussian param- 
eter [p|Jl5[ 2j| to determine the length of the time average 
as follows. The non-Gaussian parameter is defined as 



a 2 {t) 



5<r 2 (i)> 2 



- 1 



(3) 



where (•} corresponds not to the canonical average since 
the system is out of equilibrium. Instead we mean by (•} 
here and in the following an average over particles and 
initial configurations, i.e. 



n(0) 



(4) 



where the sum goes over either all A particles (to obtain 
&2a) or a ^ B particles (to obtain «2b)- This parameter 
vanishes if the van Hove correlation function 



is equal to a Gaussian u4 



(5) 



G s (r,t) 



27r(r 2 (t)) 



3/2 



cxp(-3r 2 /(2(r 2 (i)))) . (6) 



We expect Eq. (g), and therefore a 2 — 0, to be a good 
approximation for t — * (because this corresponds to the 
ballistic regime where r(t) oc v ■ t which is Maxwell dis- 
tributed p^j) as well as for t — > oo (diffusive behavior). 
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FIG. 1. Non-Gaussian parameter a^ii) (see Eq. (B)) for A 
particles. For clarity only a subset of all simulated tempera- 
tures is shown. 
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For intermediate times ai(t) ^ (see Fig. [I]). Although 
the non monotonous temperature dependence of a.% (t) in 
Fig. [l] indicates that our statistics is not very good, we use 
the time i max where 0.2(f) reaches its maximum as an es- 
timate for the characteristic time scale of mobility. (Note 
that for T < 0.25 we use i max = 10 5 .) In all following 
results which involve time averages, we use < max , if not 
otherwise stated, as the time length over which we aver- 
age. Note that i max is much larger than the microscopic 
oscillation time which is of the order of 1.0. We obtain 
thus for each temperature and for either all A or all B 
particles the distribution of df defined in Eq.(||). Fig. || 
shows for the A particles that with increasing tempera- 
ture the distribution shifts to the right and develops a 
longer tail. (Similar results are obtained for the B parti- 
cles.) The curves are zero for df < 0.0025, which reflects 
the fact that all particles are oscillating somewhat. The 
tail of P(df) extends for high temperatures to values of 
df that are twice as large as the df at the peak position, 
which shows that the dynamics is rather heterogeneous. 
We also mention that the P(df) for different T can be col- 
lapsed onto a single curve by rescaling the distribution to 
P{d 2 jm (see gj). 

Fig. y summarizes the average values of df for the A 
and B particles (df) as a function of temperature. For a 
harmonic system one would expect a linear dependence of 
(df) through the origin and over the whole range of tem- 
peratures. Since the B particles are smaller, they have a 
larger amplitude of oscillation than the A particles. For 
very small temperatures (df) increases linearly and devi- 
ates for A and B particles from a line for larger temper- 
atures. The decrease of (df ) of B particles for increasing- 
temperature at high temperatures is due to our time aver- 
age with i max . If one averages instead over the complete 
long simulation run, (df) increases monotonically and even 



FIG. 2. The distribution P(df) for the A particles. 
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FIG. 3. Temperature dependence of (df) for the A and B 
particles. The straight lines are fit to the data in the harmonic 
regime. 
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FIG. 4. (r (t)) for all A particles at the temperatures 
0.15/0.25/0.35/0.4/0.43 (solid lines) and for comparison the 
equilibrium data at T = 0.446 (bold solid line). Included are 
also the (r 2 (t)) for the fastest 5% A particles (bold dot-dashed 
lines) and a bold dashed line of slope 1. 



more than linearly. The discrepancies from a straight line 
through the origin for A and B particles at temperatures 
T > 0.25 show that anharmonic effects become important 
already at small temperatures. As suggested in [p6[ , this 
onset of anharmonicity may be related to the calorimet- 
ric glass transition. The log-log plot of the mean square 
displacement 
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FIG. 5. (r 2 (t)) for the slowest 5% A particles and B parti- 
cles. 
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n(0)r 
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where we average over all A particles (see Fig. [|), shows 
that the slopes at large times are smaller than one. There- 
fore most particles never reach the diffusive region (m = 1) 
and are trapped in their cages during the whole simula- 
tion run at least for T < 0.35. We find the same for the B 
particles where m < 0.93. If we average over only the 5% 
particles with the largest (r 2 (i cn( j = 10 5 )), the late time 
slopes are at T = 0.43 m w 1.4 for the A and m ~ 1.2 for 
the B particles. This transient behavior of m > 1 might 
be due to jump processes. Jumps are clearly visible for the 
B-particles in Fig. ||, which shows the 5% particles with 
the smallest (r 2 (t en d))- At low temperatures the slowest A 
particles are trapped at their site as can be seen in Fig. ||, 
since (r 2 (t)) < 10~ 2 over the whole simulation run. Note 
also that the slowest B particles are faster than the aver- 
age A particles. Fast B particles at T = 0.43 are reaching 
values of even (r 2 (t on d)) > 10 (see Fig. ||). 
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FIG. 6. {r 2 (t)) for the fastest 5% A particles and B particles. 



IV. SPATIAL DISTRIBUTION 



With the definition of mobile and immobile particles 
as given in the last section, we study now how they are 
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FIG. 7. Snapshot of the mobile A (white large spheres) and 
B particles (light grey small spheres) and the immobile A (dark 
grey large spheres) and B particles (black small spheres) at 
T = 0.15 and at the beginning of the production run. The 
radii were chosen for clarity and do not reflect the parameters 
of the potential. 




FIG. 8. Same as Fig. | but for T = 0.43 



spatially distributed. Fig. [7| and g show the spatial distri- 
bution of the mobile (light spheres) and immobile (dark 
spheres) at temperatures T = 0.15 and T = 0.43, and at 
the beginning of the time interval for which their mobilities 
are determined. For clarity all other 900 particles are not 
shown. In these snapshots, and similarly for all other tem- 
peratures and times, the particles are clearly distributed 
in a heterogeneous way. We therefore find dynamic het- 
erogeneity for all investigated states in the glass phase. 
The number of particles in the largest cluster p7| is at 
all investigated temperatures for mobile particles about 
30 particles and for immobile particles about 22 particles. 
These clusters are smaller than the clusters of Weeks et 
al. 0. The likely reason for this discrepancy is that we 
study a smaller system and with less good statistics. 
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FIG. 9. 3mAmA/ffAA( r ) at different temperatures. The hor- 
izontal dashed line at <? m AmA/ SAA = 1 is for the guidance of 
the eye. The inset is a comparison of <7 m AmA/<7AA (solid line) 
and <7aa(>") • 2.5 (dashed line) at T — 0.2. 

To quantify the spatial heterogeneity we plot similar to 
[ fl5[ the ratio <7 m Am A / 5a A between the radial pair distri- 
bution p8p^ ] of solely mobile particles and that of all 
particles (Fig. ^|) with 

= ^ ft -*s-d> • (8) 

In the case of randomly selected 5% particles from all A 
particles, this ratio would be one. We find however that 
this ratio is not a constant with respect to r (and similarly 
for the corresponding ratios of AB and BB), which con- 
firms the dynamic heterogeneity. Since g m AmA/ffAA > 1.0 
for distances r < 3.2, mobile particles tend to be near 
each other. We can conclude from the position of the first 
peak of 3 m Am A / 3aa (see inset of Fig. ^) that separation 
distances which are for average particles very unlikely, as 
at the left wing of <?aa, occur for mobile particles more 
often. 

We can draw similar conclusions for the immobile par- 
ticles. Fig. |l(^ shows that the ratio of the radial pair 
distribution of immobile particles to that of all particles, 
3iAiA/<?AA, is also larger than one for small distances. The 
inset, which includes $aa for comparison, reflects that also 
for immobile particles very small separation distances are 
more likely. 

V. SURROUNDING 

In the last section we found that the mobile/immobile 
particles form clusters and that the typical distances be- 
tween the particles are different from those in the bulk. 
We now address the question of the reason for mobility by 




FIG. 10. fliAiA/flAAM- The inset shows the comparison 
with <7aa- The symbols are the same as in Fig. ^. 
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FIG. 11. Number of neighbors (counting both A and B par- 
ticles) of a mobile A particle, (z m AA+B}i m comparison to the 
number of neighbors of an average A particle, (zaa+b)- 



comparing the surrounding of mobile/immobile particles 
to the one of average particles. 



A. Coordination Numbers 

To probe the immediate neighborhood of the mobile and 
immobile particles we count their number of nearest neigh- 
bors (coordination number z) where a particle j is defined 
to be a neighbor of particle i if their distance | = |fj— fj \ 
is smaller than the position of the first minimum r m j n of 
the corresponding (average) radial pair distribution func- 
tion (r min = 1.4 for AA, 1.2 for AB and 1.07 for BB, 
independent of temperatures). Fig. [ll] shows that a mo- 
bile A particle is on average surrounded by fewer particles, 
(zhiAA+b), than an average A particle, (zaa+b)- This sug- 
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FIG. 12. Number of neighbors (counting both A and B par- 
ticles) of a mobile B particle, (zniBA+B)i m comparison to the 
number of neighbors of an average B particle, (zba+b)- 

gests that one of the reasons for a particle to be mobile is 
that it is on average caged by fewer particles. The same 
is true for mobile B particles (see Fig. [l2]). However, we 
cannot make the stronger statement that any A particle 
with (z m AA+B) less than, say, 13.4 is mobile, because the 
distribution of coordination numbers is quite broad with 
a standard deviation of a « 2.0. Mobile A particles are 
furthermore surrounded by a lower than average percent- 
age of B neighbors (see Fig. |l^), because the latter trap 
A particles both energetically (cab > £aa) as well as geo- 
metrically (cab < oaa)- 

Similarly immobile particles have the property to have 
more neighbors than average particles (see Fig. |l4|) and of 
a higher percentage of B particles than usual (see Fig. jlB] 
and |lq) . Notice that the latter is true both for A and B 
particles (see Fig. |l5| and |l6|) due to the tighter packing 
with the smaller B particles. 

We also find that the percentage of mobile/immobile 
neighbors of a mobile/immobile particle is significantly 
larger than 5%, i.e., (z mAm A+mB) > 0.05 • (zhiAA+b) and 
(ziAiA+iB) > 0.05 • (^aa+b), which reflects once more the 
spatial heterogeneity discussed in Fig. and HO. 



B. Radial Pair Distribution Functions 

Next we use the radial pair distribution function to 
study the environment of the mobile and immobile par- 
ticles beyond the nearest neighbor shell. Fig. [Tl] shows 
the radial pair distribution function of a mobile A parti- 
cle with any B particle, <? m AB, in comparison with g(r) 
of any A and B particles, <?ab, at T — 0.15. We find 
that g m AB has smaller maxima and broader peaks than 
<7ab which corresponds, specifically for the first neighbor 
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FIG. 13. Fraction of B particles of all neighbors surrounding 
a mobile A particle (filled circle) and an average A particle 
(open triangle). 
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FIG. 15. Ratio of the number of B particles and all neighbors 
of an immobile A particle (filled circle) in comparison with the 
ratio of the number of B particles and all neighbors of any A 
particle (open triangle). 
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FIG. 14. Total number of neighbors of an immobile A par- 
ticle (filled circle) and an average A particle (open triangle) . 
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FIG. 16. Same figure as Fig. [15| but now for the neighbors 
of an immobile B particle. 
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FIG. 17. Radial pair distribution function of a mobile A 
particle and any B particle (solid line) and for comparison the 
corresponding radial distribution function of any A and B par- 
ticles (dashed line). The inset shows the similar g(r) for the 
immobile A particle with any B particle. Both graphs are for 
T = 0.15. 



shell, to an effectively wider cage around the mobile par- 
ticles 30 1 . The wider cage allows larger distances and 



thus larger d\ which corresponds by definition to mobile 
particles. We see the same effect for g m AA (for g m BB the 
statistics is not sufficient) and for all other temperatures. 

Immobile particles, in contrast, are surrounded by an 
effectively narrower cage, as can be concluded from the 
more pronounced peak of the first maximum (see inset of 
Fig. [It]). This is observed for all temperatures and all 
radial pair distribution functions characterizing the sur- 
rounding of an immobile particle. Notice that the change 
of the neighborhood is larger around an immobile particle 
than around a mobile particle, as the comparison of Fig. |l3| 
and Fig. |l5| and the comparison of Fig. [l?] and its inset 
show. This is probably due to our definition of mobility: 
since the distribution of d\ (see Fig. ||) is very asymmet- 
ric, 5% particles with the smallest d\ cover a much smaller 
range of d\ than 5% particles with the largest d\. Immo- 
bile particles are thus more distinct than mobile particles. 



VI. TIME SCALE OF MOBILE AND IMMOBILE 
PARTICLES 

In this section we get back to Eq. (^) , which is essential 
for the definition of mobile and immobile particles. We 
vary the time length over which we average. Specifically, 
we average over the simulation time for the long and short 
runs (see Sec. |l|), rather than using a 2 to determine a 
temperature-dependent time for calculating the average. 
We now investigate the influence of this averaging time on 
the results presented in the previous three sections. 
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FIG. 18. Total number of neighbors of any A particle 
(^AA+b) (open triangle) and of a mobile particle defined for 
the longer simulation run (z^aa+b) ( aar k filled circle) and the 
shorter simulation run (^niAA+B) (g re y diamond). 



Fig. [18] shows a comparison of the total number of near- 
est neighbors of a mobile A particle for the short runs, 



(z 



mAA+B/ 



and the long runs, (z 



mAA+B/ 



All coordina- 



tion numbers have been averaged over the 10 independent 
initial configurations, and no data from later times are 
included. The difference in (^ lAA+B ) and ( 

^mAA+B / ^ 

solely due to the different definition of mobile particles. As 
before (see Fig. |ll]) we find that the mobile particles are 
surrounded by fewer particles. For the long runs, however, 
this effect vanishes at higher temperatures. The likely rea- 
son for this decreasing difference is that we used in Eq. (Q) 
a time average over the entire long simulation run, while 
a fast particle might be fast only over some fraction of the 
simulation run. While the particle is fast, its environment 
is different than that of a regular particle. However, when 
the time average includes also times when the particle is 
not fast, then we dilute the average with environments 
that are not special. This mixing with the environment 
of average particles happens more readily at temperatures 
T > 0.4 when we approach the glass transition because the 
typical time scale over which a particle is fast is shorter 
than at lower temperatures. 

Other coordination numbers show the same behavior 



with the exception of immobile A particles (see Fig. |19| ). 
The latter distinguish themselves from the average particle 
even at high temperature and thus are immobile over the 
whole simulation run (consistent with Fig. ||). 

Another quantity which is also strongly dependent on 
the time length of the run is the distribution of 1 /df (see 
Fig. ^O] for long runs and its inset for short runs). In 
the case of long runs a double peak structure develops for 
increasing temperature. We tentatively associate the par- 
ticles in the peak with large l/df with localized particles 
and those in the peak with small 1/df with mobile parti- 
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FIG. 19. Total number ol neighbors ol an average A particle 
(zaa+b) (open triangle) and of an immobile particle defined 
for the longer simulation run (z'aa+b) (dark filled circle) and 
the shorter simulation run (^jaa+b) (s re y diamond). 
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short runs at various temperatures. 



cles. Further work is required to check this hypothesis. 

The peak at small l/df starts to dominate with increas- 
ing temperatures, which reflects that at high temperatures 
most particles have become fast at some time during the 
simulation run. 

As the figures of Sees. [V and |y| demonstrate, these 
effects of averaging over fast and regular particles can 
be avoided via appropriate choice of averaging time i max 
which we chose to be when a2(t) reaches its maximum. 
Therefore t max gives us a rough estimate about the life- 
time of fast particles. 



VII. SUMMARY AND OUTLOOK 

We investigate the mobile and immobile particles of a 
glass. Note that our definition of mobile and immobile 
is different than the definition given in jy| since we have 
the picture of a solid in mind (similar to G. Johnson et al. 
who study solid- like particle clusters [[Toll). 

We find below the glass transition temperature a clear 
dynamic heterogeneity which has previously been seen in 
simulations above T g and experiments below T g .To address 
the question why certain particles are more/less mobile 
than others, we study their surrounding. As one might 
have expected, the mobile/immobile particles are sur- 
rounded by fewer/more neighbors, forming a cage which is 
effectively wider/narrower than the one of a regular par- 
ticle. In addition mobile/immobile particles are trapped 
by fewer/more B particles, which are smaller than the A 
particles and therefore allow closer packing |31|. We ex- 
pect that, similarly, a surrounding specific to the mobility 
of the central particle might be found in the future in 
experiments. Both the dynamic heterogeneity as well as 
the particular surrounding of mobile and immobile parti- 
cles are consistent with collective behavior as it has been 
found above the glass transition temperature T g . A more 
detailed analysis below T g is left for future work. 

The characteristics of mobility show a time dependence 
which is well estimated with the time < max of the maximum 
of oii- We conclude that mobile particles are "fast" only 
for a certain time window of the simulation run whereas 
the immobile A particles seem to stay mostly immobile 
over the range of our simulation runs. This raises the 
question of a more precise criterion for the time scale of 
fast and slow processes, which we leave for future work. 
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